
cd('.')
addpath('matlab functions')

D0 = readtable('NC Data/SCC_4to5.csv');
D0.intercept(:) = 1;

OUT = struct();

for mdl = {'NOmatching','matching'}

for sbj = {'math','reading'}

    vamdl = load(sprintf('output/main/vam_%s_%s.mat',sbj{1},mdl{1}));
    vadist = vamdl.VAdistEstimates;
    controls = vamdl.controls;
    tchcoef = vamdl.tchcoef;
    outcome = vamdl.modspecs.outcome;

    gm = gmdistribution(vadist.mean',vadist.cov,vadist.pr);

    D = D0(~isnan(D0.(outcome)),:);

    numstud = 75;
    numsim = 1000;
    x = zeros([numsim numstud]);
    for n = 1:numstud
        fprintf('%d \n',n)
        for s = 1:numsim
            Ds = datasample(D,n);
            Ds.teacherid(:) = 1;
            Ds.sid = (1:n)';
            Ds.(outcome) = Ds{:,controls.Properties.RowNames}*controls.Coefficient ...
                + Ds{:,tchcoef}*random(gm,1)' + randn([n 1])*sqrt(vadist.residvar);
            vadata = genvadata(Ds,vamdl.modspecs,vamdl.controls);
            [~,vapost] = vaposterior(vamdl.VAdistEstimates,vadata);
            MOM = gmmom(vapost);
            zbar = mean(Ds{:,tchcoef},1);
            x(s,n) = zbar*MOM.cov*zbar';
        end
    end

    OUT.(sbj{1}).(mdl{1}) = struct('x',1:numstud,'y',mean(x,1));

    w = zeros([numel(vamdl.ValueAddedEstimates) 2]);
    for j = 1:numel(vamdl.ValueAddedEstimates)
        Ds = D{D.teacherid==vamdl.ValueAddedEstimates(j).teacherid,tchcoef};
        MOM = gmmom(vamdl.ValueAddedEstimates(j).vadist);
        zbar = mean(Ds,1);
        w(j,2) = zbar*MOM.cov*zbar';
        w(j,1) = size(Ds,1);
    end

    OUT.(sbj{1}).(mdl{1}).data = mean(w,1);

end

end

plotformat = {'interpreter', 'latex','FontWeight','normal','FontSize',10,'FontName','Times'};
axisformat = {'FontUnits','points','FontWeight','normal','FontSize',10,'FontName','Times'};

%-- Level Comparison --
FOOTER = .16;
DIM = [1 2];
OUTPOS = subplot_outpos(DIM,0,FOOTER,.02);

FIGRATIO = .62;

close all
pos = get(gcf,'Position');
pos(3) = sqrt((pos(3)*pos(4)/FIGRATIO));
pos(4) = pos(3)*FIGRATIO;
set(gcf,'Position',pos)

i = 1;
for sbj = {'math','reading'}
    [COL,ROW] = ind2sub(flip(DIM),i);
    subplot(DIM(1),DIM(2),i)
    plot(OUT.(sbj{1}).NOmatching.x,OUT.(sbj{1}).NOmatching.y,'-k')
    hold on
    plot(OUT.(sbj{1}).matching.x,OUT.(sbj{1}).matching.y,'--k')
    plot(OUT.(sbj{1}).NOmatching.data(1),OUT.(sbj{1}).NOmatching.data(2),'xk')
    plot(OUT.(sbj{1}).matching.data(1),OUT.(sbj{1}).matching.data(2),'ok')
    title([upper(sbj{1}(1)) sbj{1}(2:end)],plotformat{:});
    xlabel('Number of Students',plotformat{:});
    if i==1
        YLIM = [0 .05];
    elseif i==2
        YLIM = [0 0.025];
    end    
    set(gca,'box','off','YLim',YLIM,'XLim',[0 75],'OuterPosition',OUTPOS{ROW,COL},axisformat{:});
    i = i + 1;
end
hL = legend({'Level-Only Model','Matching Model','Data Average (Level-Only)','Data Average (Matching)'} ...
    ,'Orientation','horizontal','NumColumns',2,plotformat{:});
set(hL,'Position',[(1-hL.Position(3))/2 (FOOTER-hL.Position(4))/2 hL.Position(3) hL.Position(4)]);

set(gcf,'PaperPositionMode', 'manual','PaperUnits','inches'...
        ,'PaperPosition',[0 0 [1 FIGRATIO]*4.5],'PaperSize',[1 FIGRATIO]*4.5)
print('-dpdf','tables and figures/graphics/va_prec_level.pdf','-r2000')
close all  

%-- LEP

OUT = struct();

for sbj = {'math','reading'}

    vamdl = load(sprintf('output/main/vam_%s_matching.mat',sbj{1}));
    vadist = vamdl.VAdistEstimates;
    controls = vamdl.controls;
    tchcoef = vamdl.tchcoef;
    outcome = vamdl.modspecs.outcome;

    gm = gmdistribution(vadist.mean',vadist.cov,vadist.pr);

    D = D0(~isnan(D0.(outcome)),:);
    D_nonlep = D(D.lep==0,:);
    D_lep = D(D.lep==1,:);

    Z = D_lep{:,tchcoef};
    base = mean(sum((Z*vamdl.VAmomentEstimates.cov).*Z,2));

    numsim = 1000;
    numstud = [25 50 75];
    for k = 1:3
        n = numstud(k);
        x = zeros([numsim (n+1)]);
        for numlep = 0:n
            fprintf('%d \n',numlep)
            numnonlep = n - numlep;
            for s = 1:numsim
                Ds = [datasample(D_nonlep,numnonlep)
                      datasample(D_lep,numlep)];
                Ds.teacherid(:) = 1;
                Ds.sid = (1:n)';
                Ds.(outcome) = Ds{:,controls.Properties.RowNames}*controls.Coefficient ...
                    + Ds{:,tchcoef}*random(gm,1)' + randn([n 1])*sqrt(vadist.residvar);
                vadata = genvadata(Ds,vamdl.modspecs,vamdl.controls);
                [~,vapost] = vaposterior(vamdl.VAdistEstimates,vadata);
                MOM = gmmom(vapost);
                Z = D_lep{:,tchcoef};
                x(s,numlep+1) = mean(sum((Z*MOM.cov).*Z,2));
            end
        end

        OUT.(sbj{1}).(['n' num2str(n)]) = struct('x',100*(0:n)./n,'y',mean(x,1));
       
    end

end

FOOTER = .16;
DIM = [1 2];
OUTPOS = subplot_outpos(DIM,0,FOOTER,.02);

FIGRATIO = .62;

close all
pos = get(gcf,'Position');
pos(3) = sqrt((pos(3)*pos(4)/FIGRATIO));
pos(4) = pos(3)*FIGRATIO;
set(gcf,'Position',pos)

i = 1;
for sbj = {'math','reading'}
    [COL,ROW] = ind2sub(flip(DIM),i);
    subplot(DIM(1),DIM(2),i)
    plot(OUT.(sbj{1}).n25.x,OUT.(sbj{1}).n25.y,'-k')
    hold on
    plot(OUT.(sbj{1}).n50.x,OUT.(sbj{1}).n50.y,'--k')
    plot(OUT.(sbj{1}).n75.x,OUT.(sbj{1}).n75.y,':k')
    title([upper(sbj{1}(1)) sbj{1}(2:end)],plotformat{:});
    xlabel('Percent of Students LEP',plotformat{:});
    set(gca,'box','off','YLim',[0.01 .04],'XLim',[0 100],'OuterPosition',OUTPOS{ROW,COL},axisformat{:});
    i = i + 1;
end
hL = legend({'$n = 25$','$n = 50$','$n = 75$'} ...
    ,'Orientation','horizontal',plotformat{:});
set(hL,'Position',[(1-hL.Position(3))/2 (FOOTER-hL.Position(4))/2 hL.Position(3) hL.Position(4)]);

set(gcf,'PaperPositionMode', 'manual','PaperUnits','inches'...
        ,'PaperPosition',[0 0 [1 FIGRATIO]*4.5],'PaperSize',[1 FIGRATIO]*4.5)
print('-dpdf','tables and figures/graphics/va_prec_lep.pdf','-r2000')
close all  





